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The phase-separation kinetics of binary fluids in shear fiow is studied numerically in the framework 
of the continuum convection-diffusion equation based on a Ginzburg-Landau free energy. Simulations 
are carried out for different temperatures both in d = 2 and in d = 3. Our results confirm the 
C^. qualitative picture put forward by the large- Af limit equations studied in \^. In particular, the 

structure factor is characterized by the presence of four peaks whose relative oscillations give rise 
to a periodic modulation of the behavior of the rheological indicators and of the average domains 
^S) ■ sizes. This peculiar pattern of the structure factor corresponds to the presence of domains with two 

characteristic thicknesses whose relative abundance changes with time. 
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^ ■ I. INTRODUCTION 

O, . . 

f^ ' Binary fluids quenched below the dcniixing temperature exhibit an interesting off-equilibrium phenomenology that 

^^ . is nowadays reasonably well understood 0. Tipically, after an early stage when domains of the two phases are 
-4-^ formed, the kinetics proceeds by coarsening of these structures maintaining their morphology invariant. For viscous 
3 . fluids, such as some polymer blends, hydrodynamic effects can be neglected in a preasymptotic time domain and a 
^ ■ description of the dynamics in terms of diffusive processes alone is appropriate. This is the situation we focus on in 
this paper. In this case, the typical domains size R{t) grows according to the Lifshitz-Slyozov law R{t) ~ t" ~ t^/^ 
PI . This whole phenomenology is known to be correctly described in terms of the time dependent Ginzburg-Landau 
Q [ model [g[. This continuum approach, beside being often more efficient for numerical investigations, is better suited 
O ' for analytical purposes. In fact, despite an exact solution is presently lacking, successful approximate theories have 
been developed. Among these a prominent role is played by the large N expansion. Indeed, when the model is 
generalized to a vectorial order parameter with an arbitrary number TV of components, the limit N ^ 00 turns out to 
be analytically tractable and perturbation theories around it can be applied. This special limit has proven to provide 
'^ , qualitatively correct informations about the kinetic evolution contributing, in some cases, to clarify the nature of the 
■ - ' process |j]. 

When segregating systems are driven mechanically, by applying an external field or making the liquid fiow, the 
comprehension we have of the kinetics is much poorer. This is at odd with the wide technological interest of these 
systems in many application areas |^. In the case of an applied shear flow considered here, it is known that the 
evolution is profoundly changed with respect to that of a static liquid M . The most noticeable effect is the alignment 
of the domains of the two species along the fiow direction, chosen in the following as the x axis. This produces an 
anisotropy in the system and, in principle, one has to consider the typical domains sizes Rx,Ry and Rz along the 
coordinate axis separately. Modified growth exponents for R^ , Ry and Rz with respect to the case of an immobile 
fiuid are expected and a difference ax — cny = I between the power-law exponents of Rx and Ry has been reported 
|R,B|. In experiments, depending on the system considered, it has also been observed that, after an initial growth, 
the shear can stabilize the system into a stationary state characterized by a very large Rx with Ry finite, namely the 
phase separation process is interrupted p-p^ . 

Stretching of domains requires work against surface tension. Macroscopically the viscoelastic response of the fiuid 
shows up through a non- vanishing stress tensor |13jl4]. On microscopic scales it has been recognized [15l that the 
strain exerted by the fiuid can produce different effects and in particular the breakup and the recombination of 
domains. The fragmentation process is accompanied by the formation of small bubbles and more isotropic clusters. 
Then these structures, which grow by diffusion and join each other favoured by the fiow, are successively stretched 



^ 

a 



>< 



again and broken in a cyclical way. In an extended system, with many such domains, if the break-up events occur 
incoherently along the fluid we do not expect any appreciable effect, apart possibly from a global slowing of the growth 
process. Differently, if the break-up of the domains with the successive more isotropic growth is a coherent process in 
the network occurring when the stress exceeds a sort of threshold reached syncronously in all the system, one could 
expect a periodicity superimposed on the net behavior of the typical domains size and of other observables. 

In a previous paper ||l[ we have shown that a pattern similar to the latter case outlined above is indeed exhibited 
by the large- A^ model: Physical observables such as the domains radii or the rheological indicators behave for long 
times as power-laws decorated by log-time periodic oscillations. In that framework it was possible to identify this 
oscillatory pattern with a corresponding behavior of the four peaks observed in the structure factor. Moreover in the 
large- A'^ limit it was possible to show the existence of a dynamical scaling regime even in the presence of shear with 
the difference between the power law growth exponents in the flow and in the transverse directions equal to 1. This 
last observation was confirmed by a successive exact asymptotic solution of the large- TV model ||T^ . In the t ^ oo 
limit, the only accessible in this solution, however, the periodic oscillations do not survive. Then it must believed 
that these are long lasting preasyrnptotic effects. In experiments, on the other hand, a four peaked structure factor 
has been occasionally reported (ijl while the so called double overshoot of the excess viscosity observed sometimes 
[ p8| can be interpreted as being due to the oscillations of this quantity. 

Given this scenario the question of the accuracy of the large- iV limit for the description of real systems, where 
A'' = 1, must be addressed. In fact, the absence of topological defects for iV > d makes concepts like domains or 
interfaces inappropriate and the reliability of the approximation, even qualitatively, is not for granted. Concerning 
this issue important topics to be addressed are, in our opinion, i) the existence of the oscillatory pattern together with 
its counterpart of a four-fold peaked structure factor and ii) the existence of a scaling regime with specific growth 
exponents. This we do in this paper by a complete numerical investigation of the time dependent Ginzburg-Landau 
model for sheared binary fluids. A preliminary account of our results at T = in d = 2 was already published in a 
previous Letter |l3] . There the existence of structure factors with four peaks was confirmed and interpreted in terms 
of domains distributed in the system with two typical thicknesses. Here we overview and detail further these results. 
Moreover we study the effects of thermal fluctuations and the role of spatial dimensionality by comparing results for 
two and three dimensional systems. 

Regarding issue i) our results are confirmative of the large- iV prevision: four peaks of the structure factor are resolved 
and an initial oscillation of the typical domains radii and other quantities is observed. Domains with different thickness 
and structure factors with four peaks are observed also in quenching at finite temperatures and in three-dimensional 
systems. With respect to point ii) despite Rx keeps growing, CPU limitations do not allow a determination of the the 
behavior of Ry and Rz on sufficiently long timescales to detect a possible power-law behavior and the actual value 
of the exponents. Our data only provide a rough semi-quantitative agreement with the expected scaling laws which, 
therefore, cannot be inferred presently from numerical simulations alone. 

This paper is organized as follows. In Sec. (O) we introduce the time dependent Ginzburg-Landau model used 



in the simulations. In Sec. (Ill) and (IV) we present the results of the numerical simulations in rf = 2 and rf = 3 



respectively. Finally we summarize, discuss some open problem and draw our conclusions. 

II. THE MODEL 

We consider the Ginzburg-Landau free-energy 

^M = |dM|^^ + ^^'+f iv^in (1) 

where f is the order parameter representing the concentration difference between the two components. The polynomial 
terms in the free-energy density have a single-well structure when a > 0, 6 > and describe the disordered state of 
the mixture with {(p = 0). In the ordered state, for a < 0, 6 > 0, two symmetric minima are located ai ip = y^—a/b. 
These are the equilibrium values of the order parameter in the small temperature limit. The gradient-squared term 
in M) with K > takes into account the energy cost of interfaces between domains of different composition. The 
kinetics is described by the convection-diffusion equation 

^+V.(^«)^rv2^+ry (2) 

where the order parameter is coupled to an external velocity field v. Here F is a mobility coefficient and v will be 
assumed to be a plane shear flow profile, namely 



^ = "tycx (3) 

7 and ex being a spatially uniform shear rate [g|, and a unit vector in the x direction, respectively. The thermal noise 
77, which describes thermal fluctuations [p|, has zero mean and satisfies the fluctuation-dissipation relationship 

(?7(f, t)r^{f', t')) = -2Tr\/^5{f- r')5{t - t') (4) 

where T is the temperature of the fluid and (...) denotes the ensemble average. The validity of the present approach 
is restricted to systems where hydrodynamic effects can be neglected. For weakly sheared polymer blends with large 
polymerization index and similar mechanical properties of the two species, however, the present model is expected to 
be satisfactory in a preasymptotic time domain when velocity fluctuations are small py2C|] . 

Equation (0) can be cast in a dimcnsionless form after a redefinition of time, space and field scales pl|. Then we 
have chosen F = |a| = 6 = k = 1. 

The main observable for the study of the growth kinetics is the structure factor C{k,t) = {(p{k,t)(p{—k,t)), namely 
the Fourier transform of the real-space equal time correlation function. From the knowledge of the structure factor 
one computes the average size of domains in the different directions. These quantities can be defined as 

fdkC{k,t) 
j dk\kx\C{k,t) 

and analogously for the other directions. Of experimental interest are also the rheological indicators, among which 
the excess viscosity defined as PJ2^ 

^ f dk 

^V = ~-j J^kxkyC{k,t). (6) 

We have simulated equation (0) in d = 2 and d = 3 by a first-order Euler discretization scheme. Periodic boundary 
conditions have been implemented in the x and z directions (for the d — S case); Lees-Edwards boundary conditions 
^3| were used in the y direction. These boundary conditions, originally developed for molecular dynamic simulation 
of fluids in shear, require the identification of a point at {x, 0, z) with the one located at (x -I- 'jLAt, L, z), where L 
is the size of the lattice (the same in all the directions) and At is the time discretization interval. The system was 
initialized in a high temperature disordered state and the evolution was studied with a < 0. Simulations were run 
using lattices of size L = 1024 in d = 2 and L — 256 in d = 3 with Ax = 0.5, 1. The width of interfaces is given by 
\/2k |24|. We do not observe significative differences between these two choices of Ax. The results here shown were 
obtained with 7 = 0.0488, Ax = 1, At = O.Of , {(p) — 0. The two dimensional system has been studied considering 
temperatures in the range < T < 5, while only the T = case has been studied in d = 3. 

III. NUMERICAL SIMULATIONS: D=2 
A. T=0 

In this section we discuss the results of the numerical simulations of Eq. (0) without the thermal noise r]. A sequence 
of configurations of the order parameter are shown in Fig. (yj) at different values of the strain 7t. After the usual early 
stage, when domains are forming from the mixed initial state, a bicontinuous structure is observed. The distortions 
produced by the flow appear evident from 7t ~ 1 onwards. At 7t = 11 striped domains aligned with a tilt angle 9 with 
respect to the flow are observed. As time elapses 9 decreases and the structures align with the flow. In the meanwhile 
nonuniformities are formed in the system: Regions with domains of different thickness can be clearly observed. Small 
bubbles are also present, originated from the fragmentation of strained domains. 

A more accurate analysis of the spatial properties can be achieved by means of the structure factor, shown in 
Fig. (ph. Initially C{k,t) exhibits an almost circular shape, corresponding to the early stage. Then the presence of 
shear deforms C{k,t) into an ellipse and modifies the profile of the edge until, for 7t > 1, four peaks can be clearly 
observed. The presence of a peak in the structure factor is generally interpreted as the signature of a characteristic 
length in the fiuid. Here, due to the anisotropy, for each peak one associates one length for each spatial dimension. 
Since the peaks are related by the symmetry fc — > — fc one concludes that there are two distinct characteristic lengths 
for each space direction. This corresponds to the observation of thinner and thicker domains in Fig. (|l|). The relative 
height of the peaks in one of the two specular foils in which C{k, t) results to be separated starting from 7t ~ 4, is 



better seen in Fig. (Q). For a better view of C(fc, t) we have enlarged differently the scales on the k^ and ky axes. Here 
one observes that the two peaks dominate alternatively at the times 7f = 11 and "ft = 20. In a situation like that at 
jt = 11, the peak with the larger ky dominates, indicating that the stretched thin domains are more abundant. When 
the strain is further increased their rupture makes the contribution of the thick regions more important; this causes 
the other peak of C{k,t) to dominate. This mechanism is very reminiscent of what happens in the large- A'^ model, 
where the recurrent prevalence of the peaks is shown to reproduce periodically in time up to much longer times then 
those reachable in the simulations here presented. 

The mechanism of stretching and break-up which characterizes the domains evolution produces an oscillatory 
pattern in the typical size of ordered regions Rx and Ry. This is shown in the first picture of Fig. (^) [^. Ry reaches 
a local maximum when C{k, t) is of the form of Fig. (j^) at -ft = 20. As already observed thick domains are more 
abundant at this time, as one can see in Fig. (|l]) at "ft — 20. In the large- A^ model this oscillatory behavior is periodic 
in the logarithm of time and is superimposed on an average power-law increase of the length Rx,Ry Here CPU 
limitations prevent a clear check of both these predictions. This behavior reflects itself on most physical observables, 
among which the excess viscosity. This is shown in the first picture of the panel of Fig. (|g) where A77 is plotted 
against the shear strain "ft. Starting from zero in correspondence of the isotropic initial condition, A77 shows a net 
increase up to a global maximum at 7^ ~ 8. This is because the stretching of the domains caused by the shear require 
work against surface tension. This increase, however, is not monotonous and Ary shows an oscillating pattern with 
a first pronounced local maximum at "ft ~ 1.5. For larger times A77 decreases due to the ordering process in the 
fluid. Fig. da) suggests that oscillations decorates the behavior of the excess viscosity even in this regime. Another 
oscillation is completed at "ft ~ 40. This is further supported by the analogous feature in the large-iV model. We 
mention that a double overshoot of A77 is observed in experiments too ||l^ . 

B. Role of temperature fluctuations 

We have performed numerical simulations of Eq. (|2|) for several temperatures in the range — 5. Generally speaking 
we observe a similar qualitative behavior as in the zero temperature case. In particular from the analysis of the field 
configurations presented in Fig. (o) for T = 5 one again observes the anisotropic deformation of the bicontinuous 
pattern emerging after the linear regime and the inhomogeneities in the domains thickness. Thermal fluctuations are 
responsible for the increased roughness of the interfaces with respect to T = and for the thermal excitations inside 
the ordered domains, similarly to the case without shear. 

The structure factor, shown in Fig. (M), starting from "ft = A exhibits four peaks whose relative heights prevail 
recurrently. However, the difference of the heights of the peaks on the same foil is less than in the case T = 0. At 
7t = 1 instead, one has only two broad maxima along the diagonal k^ ~ ky that is splitted into four peaks only at 
later times. The same behavior is also observed at T = but for earlier times. This suggests that thermal fluctuations 
slow the evolution of C(fc, t). 

Quantitative differences with the case T — Q, possibly related to the temperature-induced slower evolution, are 
exhibited by the characteristic sizes of the domains in the x and y directions. These quantities are compared for 
several temperatures in Fig. (H). On the basis of renormalization group |13| or scaling arguments |lq] Ry and Rx 
are expected to obey an asymptotic power law growth t" with ay = 1/3 and ax — ay -{- \. The exponent ay is the 
same as without shear and is therefore presumably caused by the evaporation-condensation mechanism responsible 
for growth in the absence of flow. On the other hand the growth law in the flow direction is induced by the shear. 
An accurate determination of these exponents is not possible on the basis of our numerical simulations. Finite-size 
effects become quite soon relevant in the phase separation process. Moreover the possible power-law behaviors are 
decorated by oscillations and it would be necessary to access much longer times in order to average over this periodic 
contribution. Due to finite size effects this, in turn, would require even larger lattices, which is too numerically 
demanding for standard workstations. Despite these limitations, by looking at Fig. (0) one can still make some 
qualitative observation. At T = one observes for 7t > 30 a regime roughly consistent with ax — ay ^ 1 even if 
ax is smaller then the expected value 4/3 and ay is not even well defined because of large oscillations. By raising 
the temperature one sees that this regime is gradually changed and the effective exponent ax is lowered until, at 
T = 5, it becomes comparable to the value 1/3 observed without shear. If one believes that the asymptotic exponents 
are not changed by the strength of thermal fluctuations, as suggested by the renormalization group argument 119] 
above mentioned and by the analogy with the case 7 = 0, then the role of the temperature is that of producing a 
preasymptotic behavior characterized by a slower growth in the direction of the flow. 

In binary mixtures without shear the temperature is known to help the segregation process increasing the value of the 
amplitude in the expression of R{t) ^ At^^^ Pipll' Here, as an effect of the increased temperature, the contribution 
to the phase separation of the usual evaporation-condensation mechanism, associated to the power exponent 1/3, 



becomes more relevant than the effects of the shear induced mechanism. This results in the observed preasymptotic 
behaviour R^ = At^'^ . At later times, even if not observed in our simulations, we expect a growth with a larger 
exponent close to 4/3. 

The above trend is reflected in the behavior of the excess viscosity shown in Fig.(g). If scaling is obeyed asymptot- 
ically, so that C{k,t) ~ R,j:Ryf{kxRx, kyRy), then, from Eq. (||) one would expect Ary to scale as the inverse domain 
volume, namely Ary ~ R~^R~^. Hence, given the above discussed behavior of i?^, Ry as the temperature is changed, 
one expects Ary to gradually cross-over from a power law t"^/^ at low temperatures to t~'^/^ at higher temperatures, 
which is rougly consistent with our data. 

Finally we make some observations about the influence of the temperature on the oscillating behavior. The role of 
the oscillations is progressively suppressed increasing T. This is observed in Fig. (1q). The most apparent effect is the 
disapparence for T > 1 of the initial double overshoot of Ary. A gradual suppression of oscillations due to thermal 
fluctuations is also observed in the plot of Rx and Ry, Fig. (Q). This is consistent with the above discussion about 
the effects of raising the temperature. 

IV. NUMERICAL SIMULATIONS: D=3 

In this section we present the results of numerical simulations of the zero temperature quench process in a three 
dimensional system. The typical structures produced in the coarsening process are shown in Fig. (^ . Domains appear 
with an interconnected cylindrical shape aligned with a tilt angle 9 that shrinks in time with the flow direction. Contour 
plots of the structure factor along the coordinate planes are presented in Fig. (PI ). On the plane kx = Q the shear 
term drops out from Eq. (^) as it can be easily checked by transforming Eq. (g) into reciprocal space. Then shear 
is ineffective along this direction and the structure factor displays the usual circular form as for a static fluid. In 
the plane ky = C{k,t) has an elliptical form with axes along kx and ky. The white spots reveal that four peaks 
with the same height are present. In the plane kz a behaviour reminiscent of the two-dimensional case is observed. 
These results can be compared and better describe the experimental measures of the structure factor of p6[ . As in 
d — 2, the periodic modulation of the relative heights of the different peaks produces an oscillatory behavior of most 
observables, as shown in Fig. (|lO|) for the characteristic sizes Rx,Ry, Rz and in Fig. ( pi] ) for the excess viscosity Ajy. 
While the oscillations in the x and y directions are out of phase, due to the recurrent prevalence of thin and thick 
domains described in the two-dimensional system, it is interesting to note that Rz has a different phase with respect 
to Ry, indicating an even reacher behavior. These results provide at least an indication of a power-law increase of Rx 
and Rz, apart from the oscillations, with an exponent larger than 1.2 for Rx- In the y direction, however, there is not 
even any evidence of growth. 

V. CONCLUSIONS 

In this article we have considered the phase separation process of a fluid binary mixture in the presence of an 
applied plane shear flow. The binary mixture is described by a continuum free-energy functional and we have studied 
numerically the corresponding time-dependent Ginzburg-Landau equation with a convective term introduced by the 
flow. Results have been presented for two-dimensional systems quenched to different final temperatures and for three- 
dimensional fluids quenched to T = 0. Simulations were carried out on lattices larger than those usually considered 
in phase separation studies in order to get the highest possible resolution for the structure factor. 

In real binary fluids the existence of a four-peaked structure factor has been reported p7| ] together with a double 
overshoot of the excess viscosity. To our knowledge, however, a coherent picture indicating a possible connection 
between these apparently different phenomena with a clear physical interpretation of their origin was not available. 
The results of this paper strongly suggest that the interplay between the four peaks of the structure factor, with 
their recurrent prevalence, gives rise to an oscillatory phenomenon which reflects itself on the main observables. On 
these bases we conjecture that the double overshoot observed in experiment can be interpreted as the first part of 
an oscillatory pattern superimposed on the global trend of the excess viscosity and that other local maxima of this 
quantity could be observed in an experiment on longer time scales. 

Our simulations suggest that the physical origin of the oscillating behavior is related to the presence of two kind 
of domains, characterized by different sizes. The kinetics proceeds by stretching the thicker and breaking the thinner 
ones , so that a cyclical prevalence of one of the two types occurs. This is reflected on the properties of the structure 
factor and other quantities. These results agree at a semi-quantitative level with the pattern observed in the large- 
N approximation |l[, indicating that the global picture put forward by the N = oo theory is reliable, despite the 
different physical character of the large- A^ system without stable topological defects. This predictiveness can be used 



to infer that the observed oscillations are periodic in the logarithm of time, and that the phase separation proceeds 
via dynamical scaling, so that the characteristic sizes of the domains grow according to power-laws whose expected 
exponents have been computed in ||l6|,|l9[. However these conjectures could not be clearly addressed by our simulations 
with the available CPU resources, because much larger system should be considered on much longer times. 
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FIGURES 

FIG. 1. Configurations of a portion of 512 x 512 sites of tfie wfiole lattice are shown at different values of the strain "ft at 
T = 0. The X axes is in the horizontal direction. 

FIG. 2. The structure factor is shown at different values of the shear strain 'yt aX T = 0. The axes k^ and ky are on the 
horizontal and vertical direction, respectively. 

FIG. 3. The same structure factor as in Fig. (Q) is shown in a three-dimensional plot at 'yt — 11 and 'yt — 20 in order to 
better illustrate the different relative heights of the peaks at different times. The axes kx and ky are on the horizontal and 
vertical direction, respectively. 

FIG. 4. The average size of domains in the flow (upper curve) and in the shear direction (lower curve) are plotted as a 
function of the shear strain "ft at different temperatures. The straight lines have slope 4/3 at T = 0, 1, 3 and 1/3 at T = 5. 

FIG. 5. The excess viscosity is plotted as a function of the shear strain "ft at different temperatures. The straight lines have 
slope -5/3 at r = 0, 1, 3 and -2/3 at T = 5. 

FIG. 6. Configurations of a portion of 256 x 256 sites of the whole lattice are shown at different values of the strain "ft at 
T = 5. The X axes is in the horizontal direction. Grey-scaling from black to white corresponds to values of ip from -1 to 1. 

FIG. 7. The structure factor is shown at different values of the shear strain 7i at T = 5. The axes k^ and ky are on the 
horizontal and vertical direction, respectively. 

FIG. 8. The interfaces between domains of different compositions of a 128'' section of the simulated system are shown at 
T = 0. The a;-axis points from left to right, the y-axis points into the foil and the z-axis from bottom to top. 

FIG. 9. The sections k^ = Q, ky = Q and fcz = of the structure factor are shown at ft = 1, 3, 1, respectively, and T = Q. 

FIG. 10. Evolution of the average domain size in the x (A), y (o) and z direction (•) at T = 0. The straight line has slope 

4/3. 

FIG. 11. The excess viscosity is plotted as function of the shear strain "ft at T = 0. The straight line has slope —5/3. 
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